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ABSTRACT 

The cold classical population of the Kuiper belt exhibits a wide variety of unique physical char- 
acteristics, which collectively suggest that its dynamical coherence has been maintained through out 
the solar system's lifetime. Simultaneously, the retention of the cold population's relatively unexcited 
orbital state has remained a mystery, especially in the context of a solar system formation model, that 
is driven by a transient period of instability, where Neptune is temporarily eccentric. Here, we show 
that the cold belt can survive the instability, and its dynamical structure can be reproduced. We 
develop a simple analytical model for secular excitation of cold KBOs and show that comparatively 
fast apsidal precession and nodal recession of Neptune, during the eccentric phase, are essential for 
preservation of an unexcited state in the cold classical region. Subsequently, we confirm our results 
with self-consistent N-body simulations. We further show that contamination of the hot classical and 
scattered populations by objects of similar nature to that of cold classicals has been instrumental in 
shaping the vast physical diversity inherent to the Kuiper belt. 



1. INTRODUCTION 

The quest to understand the origins of the solar sys- 
tem dates back centuries. The last two decades, how- 
ever, have seen a renewed interested in th e problem, as 
the discovery of the Kuiper belt (Jewitt & Luu 1993) 
has provided important new clues to the physical pro- 
cesses that took place during the early stages of our so- 
lar system's evolution. The continued acquisition of new 
information gave rise to a m ultitude of new formation 
models (see Mor bidelli et al.| (^008 ) for a comprehensive 
review). Among the newly proposed scen arios, an in- 
sta bilit y model, termed the ^^Nice" model (|Tsigan is~et 




aI] |2QQ5| [Gomes et aL]|2QQ5[ jMorbidelh et al ||2UU5| ), has 



Eeen particularly successful in reproducing the obs erved 
properties of p lanetary orbits and the Kuiper belt (Lev- 
|ison et al.| |2QQ8| ). 

Within the context of the narrative told by the Nice 
mode l, planets start out in a mult i- resonant configura - 
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tion dMorbidelh et al.||2007 

and, driven by planetesimal s cattering, begm migra^ 
ing divergently ( Fernandez fc Ip| |1984). Eventually, the 
planets encounter a low-order niean motion resonance, 
which results in a transient period of instability. Dur- 
ing this period, the ice giants scatter outwards and settle 
roughly onto t heir current semi-major axes but with high 
eccentricities (Tsiganis et al. 2005 Thommes et al. 2008). 
Neptune's excited eccentricity gives rise to a chaotic sea 
between its exterior 3:2 and 2:1 mean motion resonances 
(MMRs), allowing planetesimals to ran dom- walk into the 
"classical" region ( Levison et al.|200"8 ). Subseque ntly, as 
the plane ts circularize due to dynamical friction (Stew- 



art & Wetherill 1988), the scattered and resonant popu- 
lations of the Kuiper belt are sculpted. 

An outstanding problem within the Nice model lies 
in the formation of the cold classical population of the 
Kuiper belt, which is the central theme of this study. The 
cold population is distinctive from the rest of the Kuiper 

|kbatygin@gps.caltech.edu| 



belt in a number of ways. First and foremost, as its name 
suggests, the orbital distribution is dynamically unex- 
cited. When Neptune scatters planetesimals, it tends to 
pump up their inclinations to tens of degrees. Yet the 
cold population resides on nearly co-planar orbits, with 
inclinations not exceeding ~ 5 deg ( jBrown||2001[ [Glad-] 
man et aTlpOOS ). The eccentricities of the cold popula- 



tion, on average, also tend to be diminished in compari 
son with the hot population, but the division there is not 
as apparent. Figure 1 shows the eccentricities of the cur- 
rent aggregate of observed Kuiper Belt objects (KBOs) 
between 30 and 60AU. Cold classical objects, whose in- 
clinations are below 5 deg are plotted as black dots, while 
all other objects with inclinations above 5 deg are plotted 
as blue dots. Note that cold population's eccentricity dis- 
tribution is not monotonic in semi-major axes. Between 
42AU and 45AU, planar KBOs have roughly isentropic 
eccentricities. However, low-eccentricity objects progres- 
sively disappear beyond 45 AU. We refer to this feature 
of the Kuiper belt as the "wedge" (see Figure 1). 

A second distinction is the colors of cold classical 
KBOs. In general, the Kuiper belt exhibits a vast di- 
versity of colors, from neutral gray to deep red. Within 
this range, cold classical KBOs readily stand out as 



clump of exclus ively red material (Trujillo & Brown 2002 



Lykawka fc Mukai^^2005} . In a similar manner, the size 



distribution of the cold population differs significantly 
from that of the hot classical population f Fraser et al 



2010). Finally, the fraction of binaries present ~in th e 



cold population is uniquely large ( Stephens & Noll 2006) 
Moreover, it has been shown that the wide binaries of the 
cold population in particular, would have been disrupted 
by encounters with Neptune ( Parker fc Kavelaars||2010 ), 
and thus must have never been scattered. While it is 
difficult to interpret each of these observational facts as 
conclusive evidence for a particular history, their coher- 
ence suggests that the cold classicals are a unique pop- 
ulation whose dynamical similarity has been maintained 
through the dramatic evolution of the outer solar system 
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Fig. 1. — Semi-major axis vs. eccentricity of he observed Kuiper 
belt. The black points denote objects with inclinations below i < 
5deg i.e. the cold classical population. The blue points represent 
all other objects with i > 5deg. The filled curves represent the 
scattered disk region and the major mean-motion resonances are 
shown as solid lines. The triangle, adjacent to the 2:1 MMR depicts 
the wedge structure, inherent to the cold classical population. 



dMorbidelli fc BrQwn|[2QQ4 ) . 

A number ot tormation mechanisms for the cold classi- 
cal population have been suggested. Within the context 



of a smooth migratio n scenario (Malhotra 1995 ; Murray- 
Clay fc Chiang|2QQ5[[Hahn fc M alhotr a|2QQ5| ), a primor- 



dially cold population can m principle escape dynamical 
excitation. However, other drawbacks of the smooth mi- 
gration scenario, such as the inability to reproduce sec- 
ular architecture of the planets and difficulties in form- 



ing the hot classical belt, rende r i t unlikely (jMorbidelli 
' Brasser et aLj|2QQ9[ ). [Levison et al.l ( |2UU8[) 

similar emplacement history tor the cold 



et al.||2QQ9[ 

advocated a similar emplacement history 
classicals as the hot classicals (i.e. via MMR overlap). 
The cold population that is produced in such simulations 
however, is not cold enough and not p hysically distinct 
from t he hot population. Subsequently, |Morbidelli et al.| 
(|2008") showed that if a local, col d population is imple- 
mented into the orbital solution of 'Levison et al.|p008), 
it will have the same orbital distribution as the implanted 
population after the instability, so the problem remains. 
Thus, no coherent picture of the formation of the cold 
population exists. 

In this work, we show that in-situ formation of the cold 
population is consistent with an instability model and all 
observed dynamical properties of the population, includ- 
ing the low inclinations and the wedge (shown in Figure 
1) can be formed. The aim here is not to replicate the 
Kuiper belt and the orbits of the planets in a complex 
N-body simulation, but rather to identify the dynami- 
cal processes responsible for the sculpting of the region. 
The plan of our paper is as follows: in section 2, we con- 
struct an analytical model for secular excitation of a pri- 
mordially unexcited belt, and thus derive the conditions 
for retention of dynamically cold orbits. Moreover, we 
show that the now-fossil wedge is a result of a temporary 
slow-down in orbital precession. In section 3, we perform 
self-consistent N-body simulations that conffim our an- 
alytical results. Numerical simulations show that while 
the cold population can remain undisrupted, similar ob- 
jects immediately interior to the 3:2 MMR get scattered 
all over the Kuiper belt. We conclude and discuss our 
results in section 4. 



2. SECULAR EXCITATION OF THE COLD KUIPER-BELT 

Here, we seek to develop a simplified analytical model 
that describes the long-term interactions between Nep- 
tune and an initially dynamically cold population of 
KBOs, residing between its exterior 3:2 and 2:1 MMRs, 
during a transient phase of high eccentricity. Prior to the 
instability, the planets sit in a compact configuration on 
near-circular orbits. As long as the orbital separation be- 
tween the planets and KBOs remains large, their mutual 
interactions are extremely weak, so the KBOs maintain 
their dynamically cold orbits. Consequently this period 
is unimportant to the problem at hand. 

When planet-planet resonance crossing (or some other 
mechanism) causes the instability, the gain in semi-major 
axes and acquisition of high eccentricities and inclination 
of the planets takes place on a time scale that is consid- 
erably shorter than Neptune's apsidal precession period 
(i.e. less than a million years or so). As a result, it can 
be viewed as instantaneous within the context of a sec- 
ular approximation. Thus, in an orbit-averaged sense, it 
is as if Neptune suddenly appears at 30 AU with a high 
e and i and begins interacting with the KBOs. Since we 
seek to show that after the transient phase of high eccen- 
tricity, the KBOs can end up on dynamically cold orbits, 
we must restrict Neptune from penetrating the region 
beyond 40 AU. This places the maximum eccentricity, 
attainable by Neptune, below Cmax < (4/3 — 1) = 1/3. 
This is however a weak constraint, since an eccentric, in- 
clined Neptune can still cause large modulations in the 
eccentricities a nd inclinations of the KB Os on a secu- 
lar time-scale ( Murray fc Dermott |199 9[) . Let us now 
develop a mathematical model lor these secular interac- 
tions. 

We begin by modeling Neptune's evolution. In our 
model, we take the mass of the cold KBOs to be negligi- 
ble, so they have no effect on Neptune's orbit (this is not 
necessarily true, at all times, for other Kuiper belt pop- 
ulations). The lack of mass in the primordial cold belt 
is a requirement for our model that brings up concerns 
about its formation. We shall discuss this in some de- 
tail in section 4. Since we seek to retain the majority of 
the local population, and we know that the mass of the 
current cold classical population is much less than that 
of the Earth, this is a reasonable assumption. The other 
planets, as well as the massive component of the Kuiper 
belt, will cause apsidal and nodal precession of Neptune's 

orbit, which we write as ^ = (tu^v) and / = ^^^at^ respec- 
tively. Note that we are only accounting for the average 
precessions. We express dynamical friction as exponen- 
tial decay of e fc i with constant time-scales Tg and r^. 
These time-scales are different, and their numerical val- 
ues in N-body simulation s tend to be of order ~ 10^ 
years (Levison et al. 2008). We neglect the modulation 
of Neptune's e fc 2 by the other planets. In other words, 
we only retain the free elements. 

In terms of complex Poincare variables (x = 
eexp(ztz7), 7/ = zexp(zl])), we can formulate the first- 
order Lagrange's equations for Neptune as follows: 

dXfi Xji dyn p Vn /-I \ 

-—=tgXn ~ir=^fyn • (1) 

at Tg at Ti 

It is trivial to show that these equations admit the solu- 
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Fig. 2. — Secular excitation of a KBO at a = 45AU, as dictated by 
equation (4). In these solutions, we chose ajv = 30AU, = 0.25, 
and Te = 4Myr. The final eccentricities (given by equation 8) 
are plotted as dots. Note that a low final eccentricity requires a 
comparatively fast precession. 

tions 

Xn = elexp[{ig-l/Te)t] yn = z^exp[(z/-l/ri)t], (2) 



where and are the initial (maximum) eccentricity 
and inclination of Neptune respectively and z = 
Here, we take = 0.25 and = lOdeg, in accord with 
results of numerical simulations (Tsiganis et aL 2005; 
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tated entirely by Neptune's evolution. In the spirit of 
Laplace-Lagrange secular theory, we only retain terms 
up to second order in eccentricity and inclination in the 
disturbing function of the KBOs to ensure a decoupled, 
analytical solut ion. The resulting first -order Lagrange's 
equations read ( |Wu fc Goldreich||2002 ) 

= iBykbo^iBnyn, (3) 



dt 



dt 



where A, A^, and are constants that depend only 
on the planetary masses and semi- major axes ratios of 
Neptu ne to KBOs (e.g. Ch.7 of [Murray fc Dermott 
( |1999| )). Note that in the free precession terms, [A,B)^ 
th^presence of other planets can also be accounted for 
with ease. 

From here, let us focus only on the eccentricity evo- 
lution, since the derivation of the inclination evolution 
follows an identical procedure. Setting the initial orbital 
state vector the KBO to zero ([x, y] =0), the solution to 
the above equation reads 

e^reAn{exp[iAt] - exp[{ig - l/re)t]) 
Xkbo • (4; 

Are - QTe-l 

The controlling parameter in this solution is Neptune's 
precession, g. Four solutions, for a KBO at a = 45AU, 
with different ^'s are presented in Figure 2. A natural 
unit of g is the ^8 eigenfrequency of the Laplace-Lagrange 
secular solution for the solar system, which physically 
corresponds to Neptune's average precession ra te in the 
current solar system (see Murray fc Dermott | ([l999)). 
Incidentally, the same unit can be used for the nodal 
recession rate in the inclination solution, since quantita- 
tively ^8 ^ — /s ^ 0.65" /7/r. As can be seen in Figure 2, 



varying g leads to dramatically different results. In par- 
ticular, if low-eccentricities are to be retained, g must 
significantly exceed ^8- 

After a sufficient amount of time, when Neptune's ec- 
centricity has decayed away (i.e. t ^ Te), the second 
exponential in the numerator of equation (4) can be ne- 
glected. Such a solution represents a precessing KBO 
with a constant eccentricity. Accordingly, the time de- 
pendence of the solution only governs its angular part. 
Since we are solely interested in the final orbits of the 
KBOs, we must extract only the radial part of the solu- 
tion. Let us write the t ^ solution as an exponential 
of an arbitrary number, ^: 



exp(0 



e=in( 



e^TeAn exp[iAt] 
At — gr — I 
c-expan 



(5) 



Solving for ^, and complex-expanding the logarithm, we 
have 



Vl + r2(^-A)2 



e^TeA^expzAt 

+ ^arg( — — ). 

i-Te{A-g) 



(6) 



The argument of the logarithm in the above equation 
is the radial part of the complex solution, which cor- 
responds to the final eccentricity of the KBO, with an 
equivalent expression for the inclination: 



Jinal 
■^kbo 



•final 



P^T A 



^l + riig-AY 



(7) 



In principle, we could have arrived at the same answer 
by complex-expanding the solution and taking the square 
root of the sum of the squares of the real and imaginary 
parts, although the intermediate expressions would have 
been considerably more messy. 

The above equations can be simplified even further by 
considering their limiting regimes. If the decay time-scale 
is much longer than the beat frequency {g — A, f — B)^ we 
can Taylor expand the equations to first order in (1/r) 
around zero. The answer then becomes independent of 
r. 
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This procedure is equivalent to assuming that r^{g — 
A)'^ > 1 or (/ — B)'^ > 1 and throwing away the 1 un- 
der the square root in the denominatoiQ It is clear from 
Figure 2, where the approximate solutions are plotted as 
big dots, that quantitative agreement with the "full" so- 
lution (equation 4) is excellent, in the parameter regime 
of interest. 

Figures 3 and 4 show the secular excitation of initially 
cold KBOs eccentricities and inclinations between the 3:2 
and 2:1 MMRs with g = 0,^8, 2^8 and 3^8 as solid lines. 

^ Alternatively, if the decay time-scales are short, we are in the 



non-adiabatic regime, where the solutions become e 
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Fig. 3. — Post-excitation (final) eccentricities in the cold region of 
the Kuiper belt. Solutions with g = O^gg, 2g8 and Sgg are presented 
as solid lines. Note that in order to retain nearly-circular orbits 
g > Sgs is required. The dashed line represents a solution where 
Neptune's precession rate is not kept constant. The shaded region 
corresponds to the scattered disk. 

These solutions suggest that if one is to retain an eccen- 
tricity below e < 0.1 and inclination below i < 5deg, 
Neptune's average orbital precession and nodal recession 
rates must have exceeded ~ 3^8 during the eccentric 
phase. The enhanced precession is primarily a conse- 
quence of Uranus. When Neptune scatters, it arrives 
somewhat closer to the sun than its current orbit and 
migrates to ~ 30AU by scattering KBOs (here, we have 
implicitly omitted this effect by stating that the coeffi- 
cients A, An^ and Bn are constant). Thus, at the 
time of scattering, the semi-major axis ratio of Neptune 
to Uranus may be lower, leading to an enhanced pre- 
cession. Additionally, the mass contained in the Kuiper 
belt may also play a role in inducing secular precession 
of Neptune. 

The solution described above gives eccentricities that 
monotonically decrease with semi-major axes. However, 
as already discussed above, the observed cold popula- 
tion exhibits a somewhat different behavior, with low- 
eccentricity objects progressively disappearing in the 
vicinity of the 2:1 MMR. This dynamically unique struc- 
ture (i.e. the wedge - see Figure 1) is an essential feature 
to any proposed formation mechanism for the cold clas- 
sicals. 

A wedge-like structure cannot be reproduced by a 
sweeping 2:1 MMR in an instability-driven formation 
model. Unlike the smooth migration scenario, where 
resonant capture is possible ([He nrar d fc Lamaitre||1983 



Malhotra 1995 Ketchum et al. 2Ull[ ), when Neptune is 



eccentric, the chaotic motion, that arises from resonant 
splitting (lWisdom|1980 ) , leads to an effective randomiza- 
tion of the^ccentricities (Quillen 2006). In other words, 
the KBOs that are temporarily captured do not form a 
coherent structure such as the wedge. An alternative sce- 
nario for formation of the wedge is one where the local 
population ends at 45 AU, and the wedge i s a result of 
an exten ded scattered disk with q ^ 40 AU (Gladman et 
al. 20081). It is unlikely, however, that in the extended 
scattered disk scenario, the low inclinations of scattered 
objects could be preserved. 

Here, we propose the formation of the wedge to be a 
consequence of secular perturbations. Thus, we seek to 
modify the above secular solution such that it yields ec- 
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Fig. 4. — Post-excitation (final) inclinations in the cold region of 
the Kuiper belt. Solutions with f = O^gg^ 2gs and S^fg are presented 
as solid lines. Note that in order to retain i < 5deg in the 42 — 
45AU region / < — 3gf8 is required. The dashed line represents a 
solution where Neptune's nodal recession rate is not kept constant. 
Note that the quantitative character of the solution here is subtly 
different from the eccentricity solution (Figure 3). This is because 
B involves a Laplace coefficient of the first kind, while A involves 
that of the second kind. 



centricities that are not monotonically decreasing with 
semi-major axes in the region of interest. As already de- 
scribed above, the controlling parameter in the secular 
solution is ^. So far, we have kept g constant. How- 
ever, since Neptune scatters numerous KBOs during its 
circularization, and the orbits of other planets (partic- 
ularly Uranus) are changing as well, one would expect 
Neptune's precession to vary considerably, in a chaotic 
manner. 

It is difficult to predict the exact nature of this varia- 
tion without a detailed calculation, so here we consider 
an extreme case as a proof of concept. Namely, we set 
g = 4^8 at all times, except r < t < l.lr, where we set 
^ = 0. Note that the precession of Neptune need not 
necessarily stop. We are choosing ^ = 0, rather than a 
diminished precession rate (such as, say g = gs) merely 
for the sake of argument. An analytical solution is at- 
tainable in a similar fashion as above, by breaking up the 
integration into 3 separate time intervals. If g is not held 
constant through out Neptune's circularization, the final 
eccentricity and inclination take on a different character. 
Qualitatively, this can be understood as follows: when 
Neptune stops precessing, it starts to induce consider- 
able oscillations in eccentricities of KBOs; however, once 
the precession becomes rapid again, the modulation stops 
and the eccentricities become frozen-in. These solutions 
are plotted as dashed curves in Figures 3 and 4. The de- 
tails of the non-monotonic solution depend on when and 
for how long Neptune's precession is halted, and change 
further if the precession is merely slowed down, rather 
than stopped. Furthermore, the dashed curves in Fig- 
ures 3 and 4 shift to larger semi- major axes if the free 
precessions of the KBOs (A, B) are enhanced. While it 
is understood that these calculations do not reproduce 
the cold classical population in detail, they do show that 
primordially unexcited objects can retain cold orbits in 
face of dynamical excitation, and coherent structure can 
be formed in the context of a purely secular solution. 

3. NUMERICAL SIMULATIONS 
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Fig. 5. — Orbital evolution of planets. The system starts out 
in the (2:1 4:3 4:3) initial condition, and undergoes a brief period 
of instability when Neptune (gray) and Saturn (Red) encounter a 
mutual 3:2 MMR. At t ~ 22Myr, Neptune's precession rate tem- 
porarily slows down and sculpts the wedge. The boxes on the right 
of the plot correspond to actual semi-major axes of the giant plan- 
ets. An evolved remnant planetesimal disk of this simulation is 
presented in Figures 7 & 8. 

Having motivated in-situ formation of the cold classi- 
cal population with analytical arguments, we now turn 
to numerical N-body simulations for confirmation of the 
above results and inclusion of omitted physics (such as 
close-encounters, mean-motion resonances, and higher- 
order secular terms in the disturbing function). In this 
study, the integrations were perfo rmed usin^ the mer- 
cury6 integration software package ( Chambers|1999 ) uti- 
lizing the "hybrid" algorithm. The disk was composed of 
two components. The massive planetesimal swarm, con- 
taining 3000 particles, resided between the immediate 
stability boundary of the initial multi-resonant configu- 
ration and ~ 35 AU. This was followed by a disk of an- 
other 3000 mass-less particles that extended to 60 AU. 
Thus we are assuming that a significant density gradi- 
ent exists, in the vicinity of Neptune's final orbit, such 
that the mass in the outer disk is insufficient to drive 
Neptune's migration. However, the numbers of particles 
were chosen due to considerations of computational cost 
and are not intended to be representative of the rela- 
tive fraction of bodies in the planetesimal disk in any 
way. The initial conditions were drawn from the eight 
multi-r e sonan t states that were identified by Batygin & 
iBrownl (12010') as being compatible with an instability 
lormation model. The planetesimals were initialized on 
near-coplanar, near-circular orbits (e ~ sinz ~ 10~^). 
The self-gravity of the planetesimal swarm was neglected 
to reduce computational cost of the experiments, as 30 
p ermutations of eac h initial condition were integrated. 



Batygin fc Brown| ( |2010 ) used the presence of scat- 
tering events between an ice giant and a gas giant, fol- 
lowed by a transient phase of high eccentricity, as a proxy 
for whether successful formation of the classical Kuiper 
belt can occur. Further constraints on the initial condi- 
tions can be placed by considering the reproduction of 
the outer solar systems' secular eigenmodes. Particular 
difficulty has been found in ensuring that the amplitude 
of Jupiters mode is larger than that of the ge mode 
( Morbidelli et al.||2009 ). Having completed all of the in- 
tegrations, we checked the relative amplitudes of the 
and the ge modes in all solutions. Surprisingly, we found 
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Fig. 6. — Orbital evolution of planets. The system starts out 
in the (5:3 4:3 3:2) initial condition, and undergoes a brief period 
of instability when Saturn (Red) and Jupiter (Black) encounter a 
mutual 2:1 MMR. The boxes on the right of the plot correspond to 
actual semi-major axes of the giant planets. An evolved remnant 
planetesimal disk of this simulation is presented in Figures 9 & 10. 

that despite a transient period of instability and gas gi- 
ant/ice giant scattering, the (3:2 3:2 4:3jj and the (5:3 
4:3 4:3) initial conditions did not reproduce the secular 
architecture of the p lanets, in neither this set nor in the 
set of integrations of 'Batygin & Brown (2010). If Jupiter 
and Saturn were indeed initially locked in the 3:2 MMR, 



as hydrody namic simulatio ns suggest fMas set fc Snell- 
grove 2001; |Morbideni &: C rida 2007; Pieren s ^ Nelson 
[2008), only the (3:2 3:2 5:4) and (3:2 4:3 4:3) initial con- 
artions are left as a viable options for the starting state 
of the solar system. 

As already discussed in section 2, interactions between 
the cold outer disk and the outer-most ice-giant are 
largely independent of the starting condition, since scat- 
tering in a successful simulation always sets the planets 
onto orbits that are close to that of the current solar 
system, but with moderate eccentricities. Consequently, 
we did not restrict our analysis to any particular initial 
condition. Out of our set of 180 integrations, in 8 cases, 
primordially cold objects were able to retain unexcited 
orbits in addition to the gas-giant eigenmodes being re- 
produced correctly. Here, we focus on two representative 
integrations: one starting from the (2:1 4:3 4:3) initial 
condition (Figure 5) and another starting from the (5:3 
4:3 3:2) initial condition (Figure 6). In both cases, the 
cold classical population is produced, but the wedge is 
only formed in the simulation that starts from the (2:1 
4:3 4:3) initial condition (although it is somewhat smaller 
than its observed counterpart). Note, that the formation 
of the wedge has little to do with the initial condition - 
rather, its production is a random process. Similarly, the 
exact degree of excitation of the cold population's incli- 
nations is sensitively dependent on the details of Nep- 
tune's evolution, which is chaotic. Thus, the fact that 
the wedge is reproduced in one simulation and the degree 
of excitation of the inclinations is reproduced in another 
are unrelated results. 

^ In our notation, each pair of numbers represents an MMR in 
the multi-resonant initial condition. For example, (3:2 3:2 4:3) 
corresponds to a starting state where Jupiter &; Saturn, as well as 
Saturn Sz Uranus are in 3:2 MMRs, while Uranus Sz Neptune are 
in a 4:3 MMR. 
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Fig. 7. — Eccentricity distibution of the remnant planetesimal 
disk of the simulation that starts from the (2:1 4:3 4:3) multi- 
resonant state (see Figure 5). The pale blue dots show objects 
that originated interior to ~ 35AU, 30Myr after the beginning of 
the simulation. The dark blue dots represent objects that origi- 
nated interior to 35AU, but are stable over SOOMyr. Green dots 
represent the test-particles that originate between 40 &; 60 AU. 
Yellow triangles represent test particles that originated between 35 
and 40 AU. A wedge that is somewhat similar to the observed one 
(see Figure 1) forms in this simulation, as a result of a temporary 
slow-down in Neptune's precession rate (see Figure 11). Note that 
in this simulation, the classical Kuiper belt region lies between ~ 37 
and ~ 45AU, as Neptune's final semi-major axis is a ~ 28.3AU. 
However, the aim here is to elucidate the physical mechanisms, 
rather than reproduce the actual Kuiper belt. The shaded region 
corresponds to the scattered disk. 



A vast majority (> 90%) of the objects in the cold 
classical region (i.e. a > 42AU) are retained in our sim- 
ulations on stable orbits. On the contrary, only about a 
few thousandth of the particles in the inner disk are em- 
placed onto stable orbits in the Kuiper belt region. This 
implies that in order to self-consistently study the for- 
mation of the Kuiper belt, N ^ 3000 is needed. Unfor- 
tunately, the required resolution is not computationally 
feasible. However, the problem can still be addressed by 
the use of "tracer" simulations, an approach already uti- 
lized in_the context of Kuiper belt formation by Levison 



(2008). 



et al 

in a tracer simulation the planets and planetesimals 
are not self-consistently evolved in time. Rather, the 
evolution of the planets is pre-loaded from a master sim- 
ulation and the planetesimals, which are treated as test 
particles, are evolved subject only to gravitational inter- 
actions with the planets. At the beginning of a tracer 
simulation, the tracer disk is initialized to have the same 
distribution as the massive component of the planetesi- 
mal disk. Consequently, at all times during the integra- 
tion, the tracer particles also have an identical orbital 
distribution to that of the massive planetesimals. Each 
simulation was seeded with iOO test particles and inte- 
grated on Caltech's PANGU super-c omputer. We em- 



ployed the Bulirsch-Stoer algorithm (Press et al. i992) 
in our tracer integrations. 

We performed 200 tracer simulations for each of the 
evolutions presented in Figures 5 & 6. This amounts to 
evolving a primordial disk of ~ 26000 particles, including 
the outer belt. After the ~ 30 Myr simulations were com- 
pleted, approximately ~ 7% of the particles that origi- 
nated interior to 35AU had semi-major axes in the range 
35 AU < a < 60 AU, shown as pale blue dots in Figures 7 - 



10 



% • 



• • A 

A 



2:1 



•A 



A. 

A ^ 
/ S • 



2:1 4:3 4:3 initial condition 



A • 



•A 



30 



35 



40 



45 



50 



55 



60 



,(AU) 



Fig. 8. — Inchnation distribution of the remnant planetesimal 
disk of the simulation that starts from the (2:1 4:3 4:3) multi- 
resonant state (see Figure 5). The pale blue dots show objects 
that originated interior to ~ 35AU, 30Myr after the beginning of 
the simulation. The dark blue dots represent objects that origi- 
nated interior to 35AU, but are stable over SOOMyr. Green dots 
represent the test-particles that originate between 40 Sz 60 AU. 
Yellow triangles represent test particles that originated between 35 
and 40 AU. 

10. We further cloned the population^of tracer particles 
in the Kuiper belt region to effectively increase the num- 
ber of implanted hot classical, scattered, and resonant 
particles by another factor of 6. The resulting Kuiper 
Belt, including the test-particles that originate beyond 
35 AU, was then evolved for an additional 500 Myr to en- 
sure that all unstable particles have time eject. At the 
end of the 500Myr, only ~ 5% of the implanted objects, 
that were present at the end of the 30Myr simulations, 
ended up on stable orbits. Consequently, the cumulative 
fraction of objects that are implanted into the Kuiper 
belt from the inner disk is ~ 0.3%. The stable objects 
are shown as dark blue dots in Figures 7-10. 

Note that in our simulations, the resonant populations 
are considerably diminished in number. This is largely a 
cost of performing self-consistent simulations with plan- 
etesimals that are unrealistically massive. Every time 
Neptune scatters a KBO, its resonances jump unrealisti- 
cally far, disturbin g the resonant KBOs, leading to their 
eventual ejection (Murray- Clay & Chiang 2006). In a 
suite of customized simulations where the instability still 
occurs, but planets are analyt ically guided to the ir final 
orbits, and gravity is softened ( Levison et al. 12008 ), Nep- 
tune's MMR's end up overpopulated. 'I'his leads one to 
believe that the true parameter regime of Neptune's mi- 
gration resided somew here betw een what is presented in 
this work and that of |Levison et al. (2008) (Morbidelli, 
personal communication) . " " " 

Although both of the integrations presented here pro- 
duce a cold classical belt, it is immediately apparent that 
the wedge is only produced in the integration that starts 
from the (2:1 4:3 4:3) initial condition, although again 
the process has little do with the choice of initial con- 
dition. Furthermore, from Figure 7, it can be readily 
inferred that the production of the wedge must be a sec- 
ular effect since the structure in this simulation extends 
beyond the 2:1 MMR, i.e. the unswept region. Note 

^ At the end of the simulations, there was only statistically sig- 
nificant structure in the a, e.i distributions. The orbital angles took 
on random values during scattering. 
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Fig. 9. — Eccentricity distribution of the planetesimal disk of the 
simulation that starts from the (5:3 4:3 3:2) multi-resonant state 
(see Figure 5). The pale blue dots show objects that originated in- 
terior to ~ 35AU, 30Myr after the beginning of the simulation. The 
dark blue dots represent objects that originated interior to 35AU, 
but are stable over SOOMyr. Green dots represent the test-particles 
that originate between 40 &; 60 AU. Yellow triangles represent test 
particles that originated between 35 and 40 AU. Note that the 
wedge does not form in this simulation, because Neptune's pre- 
cession never slows down, while it is eccentric. Note that in this 
simulation, the classical Kuiper belt region lies between ~ 38 and 
~ 46AU, as Neptune's final semi-major axis is a ~ 29AU. However, 
the aim here is to elucidate the physical mechanisms, rather than 
reproduce the actual Kuiper belt. The shaded region corresponds 
to the scattered disk. 

that owing to the enhanced free precession of the KBOs 
(due to presence of a massive Kuiper belt), the wedge 
structure is shifted to the right, compared with analyti- 
cal estimates presented in the previous section. 

In the context of these integrations, we are further able 
to confirm that the formation of the wedge is due to a 
considerable slowdown in Neptune's precession. During 
circularization in the integration that starts from a (5:3 
4:3 3:2) initial condition, Neptune precession is always 
roughly g ^ 4.7^8, while it is eccentric. On the contrary, 
in the integration that starts from the (2:1 4:3 4:3) initial 
condition, Neptune's precession rate varies considerably 
(l-%8 < ^ < 3.7^8) between 23Myr and 25Myr (Fig- 
ure 11). The presence of a mechanism for the successful 
formation of the wedge, from a local population is an 
important argument for confirmation of the in-situ for- 
mation of the cold classical population in the context of 
an instability model. 

It is noteworthy that in the results of the simulation, 
the wedge appears much less coherent, at semi-major 
axes interior to the 2:1 MMR. This is a consequence of ec- 
centric resonant sweeping. Because of Neptune's consid- 
erable eccentricity, the KBO multiplet and the Neptune 
multiplet of the resonance overlap even for small KBO 
eccentricities. This allows the KBO to randomly explore 
the phase-space occupied by both sections of the reso- 
nance. However, as Neptune's eccentricity is monotoni- 
cally decreasing, so is the phase space volume occupied 
by Neptune's multiplet of the 2:1 MMR, making capture 
impossible ( Quillen||2QQ6| ). Moreover, because of differ- 
ent precession rates, the nominal location of Neptune's 
multiplet of the resonance lags (i.e smaller semi-major 
axis) that of the KBO. Thus, if a KBO exits the reso- 
nance shortly after it enters, it tends to get transported 
closer to the sun, since it enters at the KBO multiplet 
and exits at the Neptune multiplet. The change in semi- 
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Fig. 10. — Inchnation distribution of the planetesimal disk of the 
simulation that starts from the (5:3 4:3 3:2) multi-resonant state 
(see Figure 5). The pale blue dots show objects that originated in- 
terior to ~ 35AU, 30Myr after the beginning of the simulation. The 
dark blue dots represent objects that originated interior to 35AU, 
but are stable over SOOMyr. Green dots represent the test-particles 
that originate between 40 &: 60 AU. Yellow triangles represent test 
particles that originated between 35 and 40 AU. 

major axes, however is only the resonant splitting width, 
so it is rather small ( Sa < O.IAU). This randomization 
of the orbital elements causes the inner part of the wedge 
to appear less coherent in Figure 7. 

It is finally worth noting that although KBOs that be- 
come the cold classical population are able to roughly re- 
tain their primordial orbital distribution, the objects be- 
tween 35AU and 40AU inevitably get scattered by Nep- 
tune during the instability. Indeed, in both simulations 
presented here, the scattered cold classicals (shown as 
yellow triangles in Figures 7-10) join the scattered disk 
as well as the hot classical population, while some parti- 
cles get trapped in resonances temporarily, during their 
evolution. 

The fact that these lifted objects mostly get emplaced 
onto stable orbits is suggestive that the results of intru- 
sion of inclined populations by cold-classical like objects, 
that took place during the instability, should still be ob- 
servable today. In other words, the in-situ formation 
scenario for cold-classicals presented here predicts that 
that a class of objects, occupying the same unique color 
region as the cold classicals, should be present in the 
excited populations. 

4. DISCUSSION 

In this paper, we present a self-consistent dynami- 
cal model for the evolution of a primordial cold classi- 
cal population of the Kuiper belt, in the context of an 
instability-driven formation scenario for the solar system. 
We show, from simple analytical considerations, that the 
cold belt can survive the transient period of dynamical 
instability, inherent to the planets. In order for a pri- 
mordially cold population of KBOs to maintain an un- 
excited state, the average apsidal precession and nodal 
recession rates of Neptune during the transient phase of 
instability must have been considerably faster than what 
is observed in today's solar system. Simultaneously, suc- 
cessful formation of the wedge (see Figure 1) requires 
that the apsidal precession rate drops by a factor of a 
few for a short period of time. Numerical integrations 
presented in this work confirm the results of the analyti- 
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Fig. 11. — Precession of Neptune's longitude of perihelion in the 
simulation that originates from the (2:1 4:3 4:3) multi-resonant ini- 
tial condition (see Figure 5). Most of the time, Neptune's preces- 
sion rate exceeds its current value by a factor of a few. However, 
the precession rate slows down considerably at t ~ 22Myr. The 
wedge forms as a result of the highlighted slowdown in Neptune's 
apsidal precession. 

cal calculations and reveal a particular result, where the 
formed cold population and the wedge closely resemble 
their observed counter-parts. The dynamical evolution of 
cold classicals we propose here is in close agreement with 
the uniqueness of cold classical's physical characteristics. 

In situ formation of cold classicals brings to light the is- 
sue of truncation of the classical belt near the 2: 1 MMR. 
In the chaotic capture mechanism proposed by |Levison| 
et al.| ([2008), the outer edge comes about naturally, as 



the 2:1 MMR sculpts the belt. In our solution, however, 
a cold belt that extends further out is surely possible. 
Thus, we are forced to attribute the proximity of the 
edge and the 2:1 MMR to a mere coincidence. Another 
question of interest is the fate of primordially cold bi- 
naries in the 35 — 40AU region. It is likely that many 
of these binaries will get disrupted by close-encounters 
with Neptune, although the exact fraction will depend 
on the details of Neptune's evolution. Consequently, an 
in-depth analysis of the evolution of the scattered cold 
KBOs may open up an avenue towards further constrain- 
ing the orbital history of Neptune. 

Although in-situ formation of cold classicals resolves 
a pressing dynamical problem within the Nice model, it 
gives rise to a new issue that requires attention. Namely, 
the outstanding question of importance is planetesimal 
formation beyond ~ 35 AU, given the steep size distri- 
bution of the cold classicals. In other words, how is the 
formation of planetesimals, up to ~ 200 km in size ac- 
complished in such a low-density environment? 

Although the answer to this question is by no means 
trivial, one possible solution to this pr oblem lies within 
the contex t of streaming instabilities ( |Youdin fc Good- 
|man""2005'). Streaming instabilities have already been 
suggested as the dominant formation process in the cold 
classical population, as gravitational collapse has been 
shown to yield wide binaries (Nes vorny et al. 2010^ . Im- 
portantly, in the proposed picture, planetesimal torma- 
tion is a threshold process, that "turns on," only when 



gas drag accumulates a critical amount of dust in a given 
location within the solar nebula. Thus, one can in princi- 
ple envision a system where most of the dust gets carried 
inward of ^ 35AU by gas drag, but infrequently, the dust 
surface density reaches a critical value in the outer neb- 
ula, causing a few, but sizable planetesimals to be born. 
Such a scenario would likely result in a very sharply de- 
creasing surface density profile in the outer nebula, at 
the epoch of disappearance of the gas. As a result, this 
picture would imply the existence of a steep density gra- 
dient in the primordial planetesimal disk, such as the 
one we require in our model, consistent with preventing 
Neptune's extended migration. 

Another possibility for the formation process is hierar- 
chical coagulation, where planetesimal growth is accom- 
plished by collis ions among smaller objects in a quiescent 
environment (Kenyon 2002). Particularly, it has been 
suggested that if aided by turbulent concentration, hier- 
archical coagulation could yield the desired mass of the 
cold classical population ( ,Cuzzi et al.|2010| ). In fact, even 
if the original mass of the cold belt exceeded its current 
value, erosion by collisional grinding could in principle be 
invoked to reduce the overall mass. However this process 
may prove problematic in repro ducing the observed wi de 
binary fraction of the cold belt ( Nesvorny et al.|[201 1 ) . 

Whatever the formation process for the cold classicals 
is, the results presented here have considerable implica- 
tions. First and foremost, the successful retention of the 
cold-classical population in the context of an instability- 
driven model, fixes the most significant drawback of the 
Nice model. Second, our scenario suggests that the cold 
classical Kuiper belt is the only population of objects 
in the outer solar system that has not been transported 
away from its formation site. Furthermore, assuming 
that collisional grinding has played a negligible role in 
the cold population's e volution (as suggested by the ob- 



served binary fraction ( [Nesvorny et al. 2011[) ), the cold 
classical population essentially yields the surface density 
of the solar nebula at a ^ 45AU, since the majority of 
the KBOs are retained in place. This potentially makes 
the cold classicals a unique laboratory for the study of 
surface processes as well as the chemistry of the primor- 
dial solar nebula. Third, based upon the results of the 
numerical simulations, we expect that objects that are 
physically similar to the cold classicals should be scat- 
tered throughout the Kuiper belt and as a result may 
explain the spectral similarity between cold classicals and 
some objects at higher inclinations. 

In conclusion, it appears that quantitative evaluation 
of planetesimal formation beyond ~ 35AU is required to 
draw a complete picture of the of the in-situ formation 
and evolution scenario for the cold classicals. However, 
the considerable improvement of the model for early 
dynamical evolution of the Kuiper belt presented here 
supports the overall validity of the hypothesis. 
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